/**
 * @file AMs.h
 * @brief 
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-04-26
 * 
 * @copyright Copyright (c) 2021  Linhuo
 * 
 * @par 修改日志:
 * <table>
 * <tr><th>Date       <th>Version <th>Author  <th>Description
 * <tr><td>2021-04-26 <td>1.0     <td>Linhuo     <td>内容
 * </table>
 */
#ifndef _solution_AMs_
#define _solution_AMs_

#include "TimeIntegrator.h"

class AM: public TimeIntegrator
{
private:
    Eigen::VectorXd alpha;
    Eigen::VectorXd beta;

public:
    void load_scheme(std::string _name);
    void One_step_forward();

};


void AM::One_step_forward()
{
    Eigen::VectorXd u_new(6);
    Eigen::VectorXd _remain(6);
    Eigen::VectorXd u_temp_1(6);
    Eigen::VectorXd u_temp_2(6);
    u_new << 0,0,0,0,0,0;
    u_temp_2 = u.col(s-1);
    _remain = u.col(s-1);
    
    for (int j = 0; j < s; j++)
	_remain += dt * beta(s-j) * fu.col(j) ;

    Eigen::MatrixXd I;
    I.setIdentity(6,6);
    Eigen::MatrixXd JF;
    JF.resize(6,6);
    Eigen::VectorXd b1;
    b1.resize(6);
    Eigen::VectorXd dif_u;
    dif_u.resize(6);
    dif_u << 100, 0, 0, 0, 0, 0;
    for(int i = 0; i < 5 && dif_u.norm() > 1e-8 ; i++)
    {   
        u_temp_1 = u_temp_2 ;
        JF = dt * beta(0) * Jfunc(u_temp_1) -I;
        b1 = u_temp_1 - dt * beta(0) * func(u_temp_1) - _remain;
        dif_u = (JF).partialPivLu().solve(b1); 
        u_temp_2 = dif_u + u_temp_1;
    }
    
    u_new = u_temp_2;

    for (int j = 0; j < s - 1; j ++)
    {
	u.col(j) = u.col(j+1);
    fu.col(j) = fu.col(j+1);
    };

    u.col(s-1) = u_new;
    fu.col(s-1) = func(u_new);

};

//depth 2
void AM::load_scheme(std::string _name)
{
    bool is_found = false;
    std::ifstream scheme("lmm_scheme");
    while (!scheme.eof() && !is_found)
    {
	std::string tmp;
	std::getline(scheme, tmp);
	scheme >> name;
	if (name == _name)
	{
	    is_found = true;
	    scheme >> p;
        scheme >> s;
        scheme >> dt;
	    std::getline(scheme, tmp);
	    std::getline(scheme, tmp);
	    alpha.resize(s);
	    beta.resize(p);
        alpha.setZero();
        alpha(0) = -1;
        double beta_Denominator;
	    std::getline(scheme, tmp);
	    std::getline(scheme, tmp);
        scheme >> beta_Denominator;
	    for (int i = 0; i < p; i++)
	    {
		scheme >> beta(i);
        beta(i) /= (beta_Denominator*1.0);
	    }
	}
	else
	{
	    for (int i = 0; i < 6; i++)
		std::getline(scheme, tmp);
	}
    }
    if (!is_found)
    {
	std::cout << "scheme not found!" << std::endl;
	std::exit(-1);
    }
};


#else
//DO NOTHING
#endif